† Corresponding author. E-mail:
Project supported by the National Natural Science Foundation of China (Grant No. 51476043), the Major National Scientific Instruments and Equipment Development Special Foundation of China (Grant No. 51327803), and the Foundation for Innovative Research Groups of the National Natural Science Foundation of China (Grant No. 51121004).
Reconstructing the distribution of optical parameters in the participating medium based on the frequency-domain radiative transfer equation (FD-RTE) to probe the internal structure of the medium is investigated in the present work. The forward model of FD-RTE is solved via the finite volume method (FVM). The regularization term formatted by the generalized Gaussian Markov random field model is used in the objective function to overcome the ill-posed nature of the inverse problem. The multi-start conjugate gradient (MCG) method is employed to search the minimum of the objective function and increase the efficiency of convergence. A modified adjoint differentiation technique using the collimated radiative intensity is developed to calculate the gradient of the objective function with respect to the optical parameters. All simulation results show that the proposed reconstruction algorithm based on FD-RTE can obtain the accurate distributions of absorption and scattering coefficients. The reconstructed images of the scattering coefficient have less errors than those of the absorption coefficient, which indicates the former are more suitable to probing the inner structure.
The theory of inverse problems has been widely developed in the past few decades due partly to its importance of applications in many important scientific and technological fields[1–5] and the advent of the scene of large computers and reliable numerical methods. The reconstruction of optical parameters in the participating medium is a typical inverse radiative problem encountered in retrieving highly-correlated parameters and has attracted significant attention in recent years. For instance, in many fields of engineering and technology, such as nondestructive testing,[6] optical tomography (OT),[7] infrared remote sensing,[8] information processing,[9] and combustion diagnosis,[10] the required optical parameters are unknown and need to be retrieved with the help of the probed radiative signals on the boundary and corresponding inverse theory. Among these applications, the OT has become a research hotspot in the field of the inverse radiative problem due to its unique advantages and wide potential applications. Theoretically, OT is an effective probe technique to reconstruct the distribution of optical parameters by analyzing near-infrared light intensities measured on the boundary of the participating medium. Given that the distribution of optical parameters is closely related to the internal structure of the medium, the optical image reconstructed via OT can also help to explore the internal structure (including the geometric shape, size or location). Compared with the traditional medicine imaging technologies, OT is non-invasive, portable, and capable of producing real-time images of clinically relevant parameters. In the field of biomedical science, the OT technology is considered to be a safer and cheaper supplement to the traditional imaging technologies and has received considerable attention on early tumor diagnosis.[11]
The reconstruction of optical parameters requires considering the laser–tissue interactions to obtain the information contained by the measured radiative intensity. The basic equation to describe the propagation of photons during the laser–tissue interactions is the radiative transfer equation (RTE). Thus an effective forward model needs to be established to simulate the RTE in the participating medium. Theoretically speaking, the forward model applied to OT can be classified into three different techniques according to the incident laser sources: the continuous laser (steady RTE),[12–14] the pulsed laser (time-domain RTE),[15–19] and the modulated laser (frequency-domain RTE).[20–23] Among these three models, the frequency-domain model, which can avoid the technical limitations intrinsic to the use of the time-domain model and gives some additional phase information compared with the steady model, has been regarded as the most promising technique to solve the inverse problems of OT accurately and efficiently. In particular, frequency domain reconstruction can reduce cross-talk between absorption and scattering parameters.[20] In the past few decades, many research groups have conducted considerable researches of the reconstruction of optical parameters based on the frequency-domain model. In the earlier research work,[24,25] the diffusion approximation (DA) model was usually employed as an approximation of the frequency-domain RTE (FD-RTE) to simplify the solving process and reduce the computational efforts. However, the DA model has limitations in estimating light propagation close to the sources and boundary or in the medium that contains the low-scattering region. Thus, more and more interest is now gradually turning to the light transport models based on solving the complete FD-RTE. Ren et al.[20] have employed the FD-RTE in the OT and have proved that the reconstruction algorithm based on the FD-RTE can reduce the cross-talk between the absorption and scattering coefficients significantly. Kim and Charette[21] have reported a sensitivity function-based conjugate gradient method to reconstruct the optical parameters with the FD-RTE, and proved that the frequency-domain data are more informative and better than steady data in separating internal parameters. Tarvainen et al.[26] have employed the Gauss-Newton reconstruction method in OT and obtained the reconstruction results of the case, including low-scattering regions based on the FD-RTE for the first time. Balima et al. have employed the discontinuous Galerkin finite element method[27] in OT to solve the FD-RTE, which can adapt the reconstruction algorithm to the complex geometry medium. Yamamoto and Sakamoto[28] have used a Monte Carlo perturbation calculation method to solve the FD-RTE in OT.
Although the complete FD-RTE has been used successfully in reconstructing optical parameters, the calculation of the gradient of the objective function with respect to the optical parameters is still problematic because almost all reconstruction algorithms are gradient-based methods in OT. The accuracies and efficiencies of these algorithms depend significantly on their gradient calculation techniques. In the past few years efforts have been made to calculate the gradient of the objective function, and many variant gradient-based methods have been developed to improve their performances in the sense of convergence speed and searching capability.[21,27,28] For instance, Lagrangian formalism, which can obtain the gradient via calculating the Lagrangian multipliers, has been widely used.[21,27] Additionally, the perturbation Monte Carlo technique[28] and gradient tree[29] have also been used to reconstruct the optical parameters. Although these algorithms used in OT can obtain the reconstructed images of absorption and scattering coefficients, the reconstructed results cannot show accurately the locations and shapes of the inclusions in the medium. This is due primarily to the fact that the obtained gradient lacks sufficient accuracy. To overcome this disadvantage, a more accurate reconstruction algorithm needs to be developed. The adjoint differentiation (AD) technique, which is always used in time-domain OT,[30,31] can calculate the gradient of the objective function accurately and efficiently without truncation error. However, the application of conventional AD technique requires the use of the numerical structure based on the time-step scheme, which is not contained in the frequency-domain model.
The motivation of the present study is to develop an accurate and efficient reconstruction algorithm based on AD technique to reconstruct the distributions of the absorption and scattering coefficients of inhomogeneous medium by using the FD-RTE. The remainder of this research is organized as follows. First, a finite volume method (FVM) formulation of the FD-RTE is well established and used as the forward model. Then, the objective function with the regularization term of the generalized Gaussian Markov random field (GGMRF) model is given to deal with the ill-posed problem. The multi-start conjugate gradient (MCG) method, which can speed up the convergence and enhance the accuracy of the reconstructed results, is employed to solve the inverse problem. A modified AD scheme which correlates the measured intensity with the inner optical parameters by using the collimated radiation is employed to calculate the gradient of the objective function. Finally, the simulation results of the forward model and the reconstruction results of the absorption and scattering coefficients in the inhomogeneous medium are discussed.
The FD-RTE is used as a forward model to simulate the propagation of photons in the participating medium. The FD-RTE in the direction
![](cpb_25_12_120201/cpb161059ieqn1.gif)
The intensity I within the medium is composed of the collimated intensity Ic and the diffuse intensity Id. The collimated component Ic obeys the following extinction law:
The collimated radiation Ic (
The measurable quantity used for OT is the exitance on the boundary, thus
The first and most essential stage in solving the complicated inverse radiation problem is to develop an accurate and efficient forward solution method. The FVM is employed to solve Eq. (
The left side of Eq. (
Converting the volume integral into a surface integral, the first term on the right side of Eq. (
Assuming that the radiative intensities on the Ac are uniform in each Ωm and are equal to the intensity in the center of Ac, then
![](cpb_25_12_120201/cpb161059ieqn2.gif)
![](cpb_25_12_120201/cpb161059ieqn3.gif)
Representing the radiative intensity within the control volume VP by the intensity
In each angle Ωm, the cell-surface intensity
Then equation (
![](cpb_25_12_120201/cpb161059ieqn5.gif)
![](cpb_25_12_120201/cpb161059ieqn6.gif)
![](cpb_25_12_120201/cpb161059ieqn7.gif)
Theoretically speaking, the reconstruction of optical parameters is an inverse radiative problem that involves the estimation of the optical parameters inside the medium by measuring the exitance intensity on the boundary of the medium. The progress of reconstruction essentially makes the estimated exitance the measured exitance by changing the estimated optical parameters gradually. The inverse problem can be solved by finding the minimum value of the objective function.
The objective function can be defined as a normalized-squared error between the predicted and measured exitance as follows:
![](cpb_25_12_120201/cpb161059ieqn8.gif)
![](cpb_25_12_120201/cpb161059ieqn9.gif)
The reconstruction of optical parameters is equivalent to the minimization of the objective function. As a simple and effective optimization technique, the CG method is used to solve the minimization problem by iteration. During the iterations, the updating directional vectors of the optical parameters are determined by the gradient of the objective function, which can be written as
The convergence speed of the CG method usually slows down with an increase of iteration times k. In particular, when
![]() | Fig. 1. Flowchart of the MCG method.[30] |
For the application of the MCG method, the core problem is to obtain the gradient dF(
Differentiating Eq. (
For the time-domain model, the derivative
Given that the objective function F is only the explicit function of the radiative intensity at the measured positions, the relation between the F and the inner radiative intensity must be established. The core significance of Eq. (
Obviously, the derivative dF/dIi calculated by Eq. (
The derivative dSc/d
A one-dimensional uniform and isotropic slab medium (see Fig.
As shown in Fig.
![]() | Fig. 3. The comparison between FVM and DFEM for (a) the hemispherical reflectance and (b) the hemispherical transmittance. |
Next a two-dimensional uniform medium (see Fig.
First, the number of solid angles is kept constant, and the FD-RTE is simulated by the forward model with grids of 11 × 11, 21 × 21, 31 × 31, and 41 × 41. Next, the number of grids is kept constant, and the FD-RTE is simulated with different quantites of the azimuthal and polar angles, respectively. The results (see Fig.
The performance of the reconstruction algorithm is tested for retrieving the different distributions of the absorption and scattering coefficients of the participating medium with the anisotropic factor of g = 0.9. The frequency-domain measured exitance of the infrared light passing through a heterogeneous medium with inclusions are simulated with the forward FVM method. The detector readings are generated on the boundary of the medium. As shown in Fig.
To compare the overall accuracy of the reconstructed images the normalized root-mean square error (NRMSE) is introduced and defined as[39]
![](cpb_25_12_120201/cpb161059ieqn13.gif)
![](cpb_25_12_120201/cpb161059ieqn14.gif)
To test the reliability of the proposed algorithm, the inclusions in the heterogeneous medium are set to be different locations and shapes. The absorption and scattering coefficients of the background medium are assumed to be μa = 5.0 m−1 and μs = 10.0 m−1. The absorption and scattering coefficients in the high extinction region are μa = 10.0 m−1 and μs = 15.0 m−1, and the absorption and scattering coefficients in the low extinction region are μa = 1.0 m−1 and μs = 5.0 m−1. The media with two square inclusions, two oval inclusions and ring inclusions are investigated in Test 1, Test 2, and Test 3 respectively.
The convergence process of the objective function in Test 1 is given in Fig.
As shown in Figs.
![]() | Fig. 8. Real distributions of (a) absorption and (c) scattering coefficients and the reconstructed distributions of (b) absorption and (d) scattering coefficients for Test 1. |
![]() | Fig. 9. Real distributions of (a) absorption and (c) scattering coefficients and reconstructed distributions of (b) absorption and (d) scattering coefficients for Test 2. |
![]() | Fig. 10. Real distributions of (a) absorption and (c) scattering coefficients and reconstructed distributions of (b) absorption and (d) scattering coefficients for Test 3. |
![]() | Table 1. NRMSEs of the reconstructed images in three tests. . |
Measurement error is inevitable in the actual experimental measurement. To illustrate the performance of the proposed reconstruction algorithm, the random error added in the measured data is considered. The measured value with random error is obtained by adding the normal distributed error to the exact measured value:
![](cpb_25_12_120201/cpb161059ieqn15.gif)
To test the robustness of the algorithm, 3%, 5%, and 10% random errors are added in the measurement data of Test 1, separately. The values of NRMSE (see Table
![]() | Table 2. The NRMSEs of the reconstructed images with considering measurement errors. . |
An accurate reconstruction algorithm for OT based on the FD-RTE is developed in the present research. The forward model using the FVM formulation of the FD-RTE is employed to predict detector readings on the surface of a medium, with the given absorption and scattering coefficients. An objective function with a regularization term is defined by the predicted detector readings and simulated measurement detector readings. The GGMRF model is used as the regularization scheme to overcome the ill-posed nature of the inverse problem. The MCG method is employed as an optimization technique to reconstruct the absorption and scattering coefficients. A modified AD scheme using the collimated intensity is employed to calculate the gradient of the objective function with respect to optical parameters. Considering the test results, the following conclusions can be drawn.
1 | |
2 | |
3 | |
4 | |
5 | |
6 | |
7 | |
8 | |
9 | |
10 | |
11 | |
12 | |
13 | |
14 | |
15 | |
16 | |
17 | |
18 | |
19 | |
20 | |
21 | |
22 | |
23 | |
24 | |
25 | |
26 | |
27 | |
28 | |
29 | |
30 | |
31 | |
32 | |
33 | |
34 | |
35 | |
36 | |
37 | |
38 | |
39 |